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Abstract 

Chiral perturbation theory predicts that in quantum chromo dynamics (QCD), light dynamical 
quarks suppress the gauge-field topological susceptibility of the vacuum. The degree of suppression 
depends on quark multiplicity and masses. It provides a strong consistency test for fermion for- 
mulations in lattice QCD. Such tests are especially important for staggered fermion formulations 
that lack a full chiral symmetry and use the "fourth-root" procedure to achieve the desired number 
of sea quarks. Over the past few years we have measured the topological susceptibility on a large 
database of 18 gauge field ensembles, generated in the presence of 2 -|- 1 flavors of dynamical asqtad 
quarks with up and down quark masses ranging from 0.05 to 1 in units of the strange quark mass 
and lattice spacings ranging from 0.045 fm to 0.12 fm. Our study also includes three quenched 
ensembles with lattice spacings ranging from 0.06 to 0.12 fm. We construct the topological suscep- 
tibility from the integrated point-to-point correlator of the discretized topological charge density 
FF. To reduce its variance, we model the asymptotic tail of the correlator. The continuum extrap- 
olation of our results for the topological susceptibility agrees nicely at small quark mass with the 
predictions of lowest-order SU(3) chiral perturbation theory, thus lending support to the validity 
of the fourth-root procedure. 

PACS numbers: 11.15.Ha, 12.38.Gc, 12.38.Aw, 12.39.Fe 
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I. INTRODUCTION 



The rich topological structure of the QCD vacuum is known to be responsible for many 
interesting nonperturbative effects, such as the chiral anomaly and chiral symmetry breaking, 
instantons, and the large mass of the r( meson. Among the wide variety of ways of looking 
at these phenomena, one may consider the effect that topological charge has on the kernel 
of the Dirac operator. It has broad implications. For example, it is intimately connected 
with the value of the chiral condensate [1]. 

The topological susceptibility Xt characterizes the tunneling rate between topologically 
distinct vacua by instantons and shows up in low energy phenomenology through the Witten- 
Veneziano formula [2, 3] and in chiral perturbation theory. A gauge configuration with 
topological charge v requires at least v fermionic zero-modes of the Dirac operator. The 
effect of quark mass on the topological susceptibility can be seen by separating the fermion 
determinant for a particular gauge field configuration into zero and non-zero modes. For N f 
flavors we have [1, 4] 

AT/ Nf 

ndet(^ + m;) = n 

/=i /=i 

where A is the imaginary part of the eigenvalue of Jj). Thus gauge configurations of non- 
trivial topology tend to be suppressed as any one of the quark masses approaches zero. 
However, this effect is compensated at increasing volume by a growing tendency of gauge 
field fluctuations to produce nontrivial topology. Chiral perturbation theory tells us [1] that 
the outcome of the competition is controlled by the parameter x = VT,m', where S is the 
chiral condensate, V is the Euclidean space-time volume, and m' is the reduced mass 

1/m' = l/mi + l/m2 + ... . (2) 

When at least one quark mass gets small at fixed volume (the "epsilon" regime, x <^ 1), 
gauge configurations with nontrivial topological charge are strongly suppressed, as implied 
by Eq. (1). In the physical regime, in which x ^ 1, which is the case for our study, 
topologically nontrivial configurations are not suppressed. Instead, leading order chiral 
perturbation theory requires that the mean squared topological charge be equal to the 
parameter x: 

{u^} = VEm', (3) 
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where the angle brackets represent an average over gauge fields. Thus the topological sus- 
ceptibility, 

Xt = {u^)/V = Sm', (4) 

remains finite in the large-volume limit. Even so, it is still suppressed as m' — 0. 

While lattice simulations of QCD have enjoyed considerable success in recent years, with 
errors on hadronic spectroscopy computations at the 1-2% level, simulations have struggled 
to reproduce this dependence of Xt on both mj and A'^^, until recently. This progress has 
come with improvements in lattice fermion technology, which has given much more control 
over chiral symmetry and lattice artifacts. 

In this article we present results for the dependence of Xt on the quark mass (through the 
taste-singlet pion mass) using improved staggered fermions (asqtad formulation). Descrip- 
tions of the asqtad formulation have been given elsewhere [5]. To eliminate contributions 
from unwanted fermion doublers, the staggered formulation takes the fourth root of the 
fermion determinant ^ det [If) + mf\ for each quark ("fourth-root procedure"), which may 
raise questions about fiavor counting. For a discussion of the issues, please see [5] and ref- 
erences therein. The primary purpose of our study, then, is to test the ability of the fourth 
root procedure to produce the correct number of sea quarks. Since the topological sus- 
ceptibility is measured directly on the gauge field configuration without the involvement of 
valence quarks, it is directly sensitive to sea quark effects. We will show that the continuum 
extrapolation of our results agrees well with lowest-order SU(3) chiral perturbation theory. 

This article summarizes results of calculations carried out over the past few years on 
ensembles with (2 + 1) fiavors of asqtad quarks as they were being generated (see the Ap- 
pendix). We continue to use the methodology of our previous work at larger lattice spacing 
and quark mass [6-8] with some refinements which appear here. The key features of our 
approach are these: 

1. obtaining the square of the topological charge from the integral of the two-point cor- 
relator of the topological charge density. 

2. reducing the variance of the integral by modeling the asymptotic form of the correlator 
in terms of known hadronic contributions, and 

3. analyzing the quark-mass and lattice-spacing dependence of the resulting susceptibility 
in terms of predictions from rooted staggered chiral perturbation theory. 
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In the following section, we discuss the details of our method for calculating the topolog- 
ical susceptibility on the lattice. We present our results and analysis in Sec. III. Finally, we 
comment on our results in the Conclusions, Sec. IV. The Appendix lists the parameters of 
the gauge field ensembles used in this study. 

II. METHODOLOGY 

A. Definition of the topological susceptibility 

We introduced the topological susceptibility in Eq. (4) as the mean squared charge per 
unit volume: Xt = (j^^)/^- The net topological charge v is the integral over Euclidean 
space-time of the topological charge density. 



The susceptibility is then given by the integral of the correlator of the charge density, 
provided the integral is well defined. 



where r = \x\. Because the exponential decay of the correlator at large r is set by nonzero 
hadron masses, we see that the susceptibility is properly regarded as a local observable, 
i.e., it can be defined in terms of a correlator that has finite physical range. We use this 
definition of the susceptibility, coupled with a smeared lattice discretization of FF. 

In the continuum limit the integral definition above is problematic. The unregulated 
correlator C(r) is nonintegrable: it has a positive, divergent contact term (at the origin) 
and, close to the origin, a compensating negative ultraviolet singularity of order (up to 
possible logarithms) [9-12]. Cancellation is required in order to produce the expected 
finite integral Eq. (4). To circumvent this mathematical difficulty Liischer formulated a 
definition of the topological susceptibility in terms of a product of pseudoscalar and scalar 
densities of Ginsparg- Wilson quarks [13]. Since the definition requires computing all-to-all 
quark-line disconnected correlators, it is more difficult to implement, and, to our knowledge, 
it has not yet been put into practice. 

For present purposes we resort to the naive definition in Eq. (6) and trust that the lattice 
cutoff and a smoothed definition of p{x) regulate the compensating singularities enough over 



p(x) 



rpa rpa 

32^2 Ml'- 



(5) 




(6) 
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a range of reasonably small lattice spacings that we can test the expected suppression of 
the topological susceptibility. In our present scheme we fix the smoothing scale in lattice 
units as we take the lattice spacing to zero. Our numerical simulation provides a practical 
test of the limitations of such a scheme. If it fails, as the lattice spacing is decreased, we 
would expect to encounter a growing variance from contributions to the integral near the 
origin. This would not invalidate the method: The central value in the continuum limit is 
finite even if the variance is unbounded. It could, however, require an impractically large 
computational effort to achieve a desired accuracy as we make the lattice spacing smaller. 
We return to this question after presenting our results. 

There are a variety of lattice methods for obtaining the topological charge. The traditional 
"algebraic" method uses a lattice discretization of the density FF, constructed at each lattice 
site from appropriate closed loops of gauge links. To suppress ultraviolet noise at the cutoff 
scale, smoothing is required [14]. The Boulder discretization [15, 16], which we use in the 
present study, is a refinement of the traditional definition. It is fattened (smoothed) by first 
performing some number (we use three) of HYP smoothing sweeps [17] on the gauge field 
and then constructing the operator from the smoothed links. 

A more elegant method defines the topological charge density in terms of a chiral (e.g. 
overlap) Dirac operator D, as p{x) oc Tr[75D]^ ,j. [18, 19] (the trace is over color and spin), but 
using it directly in Eq. (6) is computationally expensive [20, 21]. For the overlap operator 
a more tractable method uses the Atiyah-Singer index theorem to relate the topological 
charge v to the net number of zero crossings of the low-lying eigenvalues of a Hermitian 
Dirac kernel from which the chiral operator is built [22]. This method was implemented in 
[23]. For the overlap operator, smoothing is inherent in the choice of the Dirac kernel from 
which the overlap action is built. 

Another promising method works with gauge configurations of fixed topology [24, 25]. 
In this case, at large distance the correlator of the topological charge density approaches 
a constant Xt/^ plus other known constants that depend on the fixed topological charge. 
One can also use a hadronic flavor- singlet interpolating operator with J^*^ = ^ as a 
proxy for FF. This method has been tested at one lattice spacing in the two-flavor case on 
configurations generated with the overlap action [26]. 

The Liischer definition [13], based on a chiral Dirac operator, replaces the integral of 
FF with the integral of a quark pseudoscalar density. The quark field from which that 

6 



density is constructed can have arbitrary mass, which sets the locahzation scale of the 
operator. The expectation value of that density is regulated with a suitable number of 
zero-momentum scalar-density insertions on the quark line. At large mass in the hopping 
parameter expansion, the operator can be expressed as a sum of gauge-link loops analogous 
to those in the Boulder discretization, which regulates the construction of FF through an 
extended discretization and HYP-smeared gauge-links. In the Boulder case the localization 
of the gauge paths is controlled by the number of smearing steps, whereas localization of 
the Liischer operator is controlled by the quark mass. Of course, the chiral properties of the 
underlying action in that case allows an arbitrary choice of scale. 

Whatever the definition, the resulting susceptibility is subject in general to multiplicative 
and additive corrections at nonzero lattice spacing [27]: 



An additive renormalization is not required for chiral actions that use the same operator 
in the fermion determinant and the measurement of the topological charge [28]. In our 
case an additive renormalization is expected. We assume that in the continuum limit M 
approaches one and A approaches zero. Since with our actions lattice artifacts appear at 
0{a^) (up to logarithms), we expect that the approach to these limits is as a? [6]. With the 
overlap method one can use the same Dirac operator for the Monte Carlo evolution and the 
measurement of topological charge. In this case the small instantons and dislocations that 
are not seen by the overlap operator, so not suppressed by a small quark mass, are then also 
not seen by the topological charge operator. In our case the Monte Carlo Dirac operator 
and topological charge operators are unrelated, so we might expect larger lattice artifacts. 

B. Predictions from chiral perturbation theory 

Our computed topological susceptibility is a function of the quark masses and the lattice 
spacing. As we have already recalled in Sec. I, in chiral perturbation theory the susceptibility 
Xt depends on the number of light quarks and their masses in leading order through 



where S is the chiral condensate to this order, m„, m^, and m<j are the masses of the up, 
down, and strange quarks, and the ellipsis represents contributions beyond the cutoff from 
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(7) 



Ts/Xt = i/rriu + l/nid + l/m^ + 



(8) 



higher quark masses and the axial anomaly [1]. We see that as quark masses vanish, the 
susceptibility must vanish. The rate at which it vanishes depends on the number of light 
flavors. 

For equal up and down quark masses we may use the Gell-Mann-Oakes-Renner relation, 
also from leading order chiral perturbation theory, to rewrite this expression as 

fl/{Axt) = 2/ml + 1/ml + ..., (9) 

where m^^ = 2m|- — is the squared mass of the flctitious pseudoscalar meson containing 
two nonannihilating quarks with masses equal to the strange quark, and in our normalization 
the pion decay constant /tt is approximately 130 MeV. 

With nonchiral lattice fermions, at nonzero lattice spacing one should instead use a version 
of chiral perturbation theory appropriate to the lattice fermion formulation. In this way some 
of the lattice discretization errors can be modeled. For staggered fermions using the fourth- 
root procedure, we use rooted staggered chiral perturbation theory (rS^PT) [29]. This theory 
has a taste multiplet of sixteen pions. Among them, only the taste singlet pion is sensitive 
to the anomaly and so enters the expression for the topological susceptibility at leading 
order. At tree level the continuum expression is modifled by replacing the pseudoscalar 
meson masses by their taste-singlet counterparts [30]: 

1/Xt = {A/f!){2/ml, + 1/mlj + 3/ml), (10) 

where the subscript / identifles the taste singlet, and through the term in mo, which is 
proportional to the t]' mass at lowest order, we have introduced an explicit anomaly contri- 
bution. The standard chiral perturbation theory expression corresponds to ttlq — )■ oo (and 
a — 0); introducing mo in Eq. (10) is phenomenological because mo is beyond the physi- 
cal cutoff scale of chiral perturbation theory. At inflnite quark mass we get the quenched 
topological susceptibility xtq, which suggests an alternative phenomenological form [31], 

= i^/m/ml, + l/mL,,) + 1/xt, . (11) 
C. Topological charge density operator 

As before [6], we use the topological charge operator of DeGrand, Hasenfratz, and Kovacs 
[15] optimized for SU(3) by Hasenfratz and Meter [16]. The operator is constructed from 
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closed ten-link paths of gauge matrices as follows: 

p{x) = j: c]Tr(l - t/,) + c?[Tr(l - f/,)]^ • (12) 

Specifically, the operator Ui is constructed from a product along a path from site x in the se- 
quence of directions {x, y, z, —y, —x, t, x, —t, —x, —z), summed over rotations and reflections, 
and the operator U2, from the directions (x, y, z, — x, t, —z, x, —t, —x, —y). Both paths lie in- 
side a 2^ hypercube. The coefficients are c\ = 0.07872507, 4 = 0.3173630, cj = -0.1888383, 
and C2 = 0.2854577. Hasenfratz et al. devised this operator to optimize a match with a 
geometric definition of topological charge on a "typical" set of gauge configurations. The 
operator also reproduces accurately the charge of an instanton, provided the instanton ra- 
dius is larger than the lattice spacing. The finer details of the construction of this operator 
are unimportant for our purposes, since in the end we take the continuum limit. 

We applied this operator to gauge configurations smoothed by three HYP steps [17]. 
From the point of view of the unsmoothed gauge field, this operation, in effect, enlarges the 
footprint of the topological charge density operator by a small amount. We have shown in 
[6] that the topological susceptibility on a coarse lattice (a ~ 0.12 fm) is constant within 
statistical errors of 8% for one to four HYP sweeps. 

D. Variance reduction method 

We calculate the topological susceptibility by integrating the topological charge density 
correlator in Eq. (6) over the lattice four-volume. In the left panel of Fig. 1 we show a 
typical correlator C(r). It is expressed in units of the Sommer parameter Tq ~ 0.454 fm 
[32]. As expected, it has a positive peak at the origin next to a negative minimum, and 
it rises to its asymptotic limit of zero from below as required by CP symmetry. To give a 
better visual impression of contributions to the susceptibility, in the right panel of Fig. 1 we 
multiply C(r) by the statistical weight factor w{r) that counts the number of lattice points 
that, by symmetry, have the same four-radius r, or, where the plotted value is binned, 
have the same range of four-radii. This is essentially a discretized version of r^C{r)dr. The 
irregular binning inherent in the discretized distance r produces the ragged appearance of the 
weighted values at small r. On the other hand, statistical fluctuations produce the ragged 
appearance at large r. The topological susceptibility in ro units is simply proportional to 
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FIG. 1. Left: Topological charge density correlation function C(r) vs. separation in units of tq. 
Right: Correlation function weighted by the volume measure. Errors are statistical and have not 
been corrected for autocorrelations. The red symbols (crosses) indicate the fitted points. The black 
curve shows the fit, which we use to replace the measured points for r > rc, the cut radius. (The 
lone symbol at the right bins all measurements for r/ro > 10.7). 

the sum of the weighted values. 

In Fig. 1, right, the substantial cancellation of positive and negative contributions at small 
r is more evident. We also see that the large distance contribution to the susceptibility is 
mostly noise. We have found that it is responsible for the bulk of the variance in the integral. 
This is to be expected. In a suitably large subvolume Vq of spacetime, we should be able 
to determine the topological susceptibility reasonably well by measuring fluctuations of the 
local topological charge uq. Consider putting together such volumes to create the total 
volume V. The overall topological charge u is then obtained as a random walk of local 
charges, so its variance grows with A^. We can measure the susceptibility in two ways: (1) 
average the locally determined {i^q)/Vq over the N subvolumes or (2) calculate {u'^)/V over 
the full volume. With the former method the error in the measured susceptibility decreases 
as I/a/ZV with increasing and fixed Vq, whereas with the latter method the error never 
improves. 

In our case the integral of the correlator C(r) replaces the sum over subvolumes. But 
we still need to eliminate noise from contributions at large r. To do so, several years ago 
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we introduced a variance-reduction method that fits the large r part of the correlator to 
its asymptotic form in Eq. (15) and then, for r > for a suitable cutoff Tc, replaces the 
numerical sum of the correlator with an integral over the fitted function as follows [7]: 

Xt= I C{r)+j Cfit(r), (13) 

Jr<rc Jr>rc 

where is chosen inside the fit range. In the present study we chose ^ 1.2 fm for all 
ensembles. In Fig. 1, right, we illustrate the fit to the large r part of the correlator and 
indicate Vc- We continue to use this method in the present work. 

E. Asymptotic fit model 

The topological charge density is a flavor- singlet operator with quantum numbers J^'" = 
^, so the asymptotic behavior of the correlator is governed by the rj and t]' mesons and, 
for sufficiently light sea quarks, by multipion states. That is 

C{r) = {p{x)p{0)) ^ A,S{m„ r) + A,,S{m,>,r) + ..., (14) 

where the A^s are overlap constants and S{m, r) is a scalar propagator with asymptotic form 

S{m,r) ^ [1 + 3/(8mr)] exp(— mr)/r^''^. (15) 

The three-pion continuum is the lightest multimeson state in this correlator. For our en- 
sembles the 7] meson is always lighter. Furthermore, the coupling of the topological charge 
density operator to multimeson states is Zweig-rule suppressed. Therefore, we ignore them 
in the present analysis. Since the topological charge density operator is an SU(3) flavor 
singlet, it couples to the flavor singlet component of the rj and rj' mesons. In the usual 
representation of singlet-octet mixing [33] , 

|?7) = cos 6 Ir/s) + sin 9 |?7o) 

|?7') = - sin^ |?78) cos6' |?7o) , (16) 

so 

Ari/Ar^, = tan^ 9. (17) 

Our statistics are insufficient for determining all the parameters of the fit function reliably. 
Instead, we model the masses of the r] and r]' and the ratio A^^/A^j', leaving only one fitting 
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parameter A^j', which simply sets the normahzation of the asymptotic form. We set = 
2m^g/3 + rri^/S for our measured lattice values of m^^ and m^r, and we fix m^/ = 958 MeV 
(its physical value) since we have not calculated it for these ensembles. Finally, we use a 
simple chiral model to fix the ratio of couplings Arj/Arj' or equivalently, the singlet-octet 
mixing angle as a function of quark masses. 

Our model is based on the mass matrix for the flavor-neutral taste-singlet mesons in 
lowest order SU(3) chiral perturbation theory [29] 



M 



rrin 



rrin 



\ 



(18) 



V ^si + ^IJ 

where Mm and Msi are masses of unmixed uu {dd) and ss meson states, and ttiq parame- 
terizes the anomaly. The isosinglet eigenvectors are 



where 



Vu 
Vs 



\r]) = Vu \uu) + Vd ddj + \ss) 
W) = + ^d\dd) + v'Jss) , 

Vd = l/N 

-{MIjj - Ml J + ml + Vd)/{2mlN) 
v'd = VN' 



(19) 



(20) 



MIj + ml - ^/d) / i2mlN') 



d = (Mj, - MIj.Y - 2(Mj, - Mlj,)mi + 9m^, 



and and A^' normalize the eigenvectors to 1. Since the flavor singlet state in this basis is 
just (1, 1, l)/-\/3, we obtain the ratio 



l\2 



Ari/Ar,' = tan 6 = {v^ + Vd + v.,) / {v'^ + v'd + v'^) 



(21) 



which we apply to the fit model of Eq. (14). To complete the model, we need the value 
of the anomaly parameter ml. We set it so that for physical values of M^jj and M|j- {i.e., 
values that give the physical masses and mg^ = J2m\- — m1), we get the standard 



phenomenological mixing angle 9 ^ —20 degrees [33]. At this "physical" point the mixing 
model also gives us = 493 MeV and m^/ = 953 MeV, reasonably close to their physical 
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values. Then for unphysical masses we use the lattice values of M^^ and M^j on each 
ensemble, always keeping mo fixed. This procedure assures that the 77 decouples as required 
in the SU(3) flavor limit = md = rris, and it provides a smooth interpolation between 
that limit and the physical limit. The taste-singlet masses M^j and M^j are obtained by 
adding measured or estimated taste splittings to the masses of the lightest members of the 
taste multiplet. Splittings are listed in Table V below. 

The model is applied to all the dynamical ensembles in this study, listed in the Appendix 
A. The resulting fit parameters are listed in Table I. The mixing parameter A^/A^r is shown 
to three digits. Apart from systematic errors in the model itself, in principle it inherits a 
statistical error from our determination of the taste-singlet masses, which, in turn depends 
on the error in the taste splitting. The last error, however, is less than 5%, small enough to 
have no effect on the mixing parameter to the number of digits reported. The remaining fit 
parameters do not depend on the taste-singlet masses. Consequently, statistical errors in the 
determination of the taste-singlet masses have negligible effect on results for the topological 
susceptibility. 

F. Asymptotic fit model for the quenched ensembles 

For the three quenched ensembles we use the same methodology, except that the fit model 
has only one mass. We fix it to the mass of the J^'" = ground state lattice glueball 
from Chen et al. [34], namely 2560 MeV. The parameters are listed in Table II. We chose 
Tc for the quenched ensembles to match our choice for the dynamical ensembles at the same 
lattice spacing. Since the quenched correlators die so quickly at large r, the contribution 
to the susceptibility for r > Tc is negligible, and the asymptotic model has no effect on the 
result. 

III. RESULTS 

We smooth the lattices with three HYP smoothing steps [17] and measure the topological 
charge density with the Boulder operator at each space-time point. We then construct the 
point-to-point correlator C(r) for every pair of points in the space-time volume. For r/a < 5 
we keep values for every displacement, and for larger r we bin data over small intervals in 
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10/5' 






am^ arrir,' xLw X^/df 






coarse 


6.85 


0.05/0.05 


0.000 


485 573 6 6 3 9/11 


6.83 


0.04/0.05 


0.010 


0.470 0.578 7.6 4.5/11 


6.79 


0.02/0.05 


0.095 


0.439 0.583 17.0 10.0/11 


6.76 


0.01/0.05 


0.166 


0.424 0.588 8.8 5.2/11 


6.76 


0.007/0.05 


0.194 


0.417 0.584 11.1 6.5/11 


6.76 


0.005/0.05 


0.215 


0.413 0.582 11.8 6.9/11 






fine 


7.18 


0.031/0.031 


0.000 


0.320 0.403 40.5 11.2/19 


7.11 


0.0124/0.031 


0.072 


0.292 0.415 43.6 12.1/19 


7.09 


0.0062/0.031 


0.128 


0.280 0.416 24.3 6.8/19 


7.085 0.00465/0.031 0.144 


0.277 0.416 33.7 9.3/19 


7.08 


0.0031/0.031 


0.162 


0.274 0.417 17.9 5.0/19 


7.075 0.00155/0.031 0.181 


0.271 0.416 26.3 7.3/19 






superfine 


7.48 


0.0072/0.018 


0.049 


0.186 0.291 52.4 16.4/29 


7.475 0.0054/0.018 


0.066 


0.182 0.291 29.8 9.3/27 


7.47 


0.0036/0.018 


0.087 


0.178 0.291 30.4 9.5/27 


7.465 0.0025/0.018 


0.101 


0.175 0.291 26.5 8.2/26 


7.46 


0.0018/0.018 


0.110 


0.174 0.292 35.2 11.0/26 






ultrafine 


7.81 


0.0028/0.014 


0.097 


0.136 0.216 29.1 14.6/30 



TABLE I. Parameters used in asymptotic fits to the (2 + l)-flavor topological charge density 
correlator. The raw is uncorrected for autocorrelations. The last column includes the correction 
as explained in Sec. Ill A. 

r. The resulting data is then fit to Eq. (14) over a range [r mm, ''"max]- We replace the raw 
data with the fit model for r > r^. The fit range is chosen to give an acceptable /df 
(corrected for autocorrelations) and to vary smoothly as a function of sea quark mass and 
lattice spacing. 
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10/5-^ aniG x^/df 



8.00 1.55 16.0/12 



8.40 1.11 9.8/11 



8.80 0.816 10.0/13 



TABLE II. Parameters used in asymptotic fits to the quenched topological charge density correlator 
and resulting values of x^/df- 

A. Monte Carlo time histories and autocorrelations 

To determine the confidence level of our fits and errors in the fit parameters, we must 
first analyze autocorrelations in Monte Carlo time. With our action and molecular dynamics 
algorithm, the total topological charge is moderately persistent in Monte Carlo time. In 
Fig. 2, we show the time histories for a range of lattice spacings for mud = 0.2ms ensembles. 
As we have noted, however, the topological susceptibility is a local observable. We can get 
a graphical sense of the autocorrelation affecting the susceptibility by considering the time 
history of the integral of the correlator 



In Fig. 3 we show the time history of this variable for the case r = 2ro for the same set 
of ensembles. Clearly the fluctuations in this quantity decorrelate much more rapidly than 
those of the total topological charge. 

We estimate the autocorrelation correction, i.e., the amount by which the naive (uncor- 
related) variance should be increased to compensate for autocorrelations. For this purpose 
we consider the integral of the correlator over the proposed fit range 



We block the data in Monte Carlo time and calculate the variance of the mean as a function of 
block size, extrapolating to infinite block size. The ratio of the extrapolated variance to the 
naive variance is the correction factor. We also sum the autocorrelation coefficients to obtain 
another estimate of the correction factor. These determinations fluctuate as a function 
of sea quark mass. We averaged them at fixed lattice spacing to obtain the correction 
factors shown in Table III. We should emphasize that the determination of autocorrelation 
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2000 4000 

traj 

FIG. 2. Total topological charge after three HYP sweeps as a function of simulation time for four 
lattice spacings and fixed sea quark masses with ratio m^^/jn^ = 0.2. Sections marked "a" and 
"b" come from different Markov chains. From top to bottom, a = 0.12, 0.09, 0.06, and 0.045 fm. 

corrections is notoriously difficult. To develop more confidence in these estimates, we should 
have considerably longer time histories. 

Our fits to the data take into account correlations in r as well. For all ensembles, mea- 
surements are taken every six or sometimes every five molecular dynamics time units. We 
do not bin data in Monte Carlo time before constructing the covariance matrix in r and 
minimizing the correlated [35]. Uncorrected errors are derived from a jackknife analysis. 
Thus the resulting x^, based on the naive covariance, must be reduced by the factor in 
Table III before estimating the confidence level. Furthermore, the naive single-elimination 
jackknife errors in the fit parameters must be increased by the square root of this factor. 

16 





0.06 




0.04 




0.02 




0.00 




0.06 




0.04 




0.02 




0.00 




0.04 




0.02 




0.00 




-0.02 






o 


0.04 






-i-J 


0.02 


X 




o 


0.00 




-0.02 








2000 4000 
traj 



FIG. 3. Contribution to the topological susceptibility for r < 2ro as a function of simulation time 
for the ensembles of Fig. 2. 

We use the same factor to adjust the error in the contribution from the raw data for r < Tc- 
B. Topological charge density correlator 

We expect the topological susceptibility to decrease with decreasing light sea quark mass. 
It is interesting to see how the topological charge density correlator itself varies with the 
light sea quark mass at fixed lattice spacing. In Fig. 4 we examine this dependence for a 
series of fine lattice ensembles (a ~ 0.09 fm) for which we have results for four ratios of the 
light to strange quark mass, 0.05, 0.1, 0.15, and 0.2, corresponding to the range 0.601 to 
1.074 in m^^rQ. In the upper panel any variation with light quark mass is evidently much 
smaller than the plot symbol size. In fact the short distance part of the correlator shows 
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spacing correction 
coarse 1.7 
fine 3.6 
superfine 3.2 
ultrafine 2.0 

TABLE 111. Autocorrelation correction factors for the various categories of lattice spacings in this 
study. The factor multiplies the naive variance. 

very little sea quark mass dependence. In the lower panel we enlarge the region around the 
minimum where a small variation is now apparent. In this region hght meson states begin 
to dominate the correlator of the gluonic operators. As the hght quark mass decreases, the 
minimum drops, thus giving a larger negative contribution to the integral. This effect leads 
to the suppression of the susceptibility. According to the model, the correlator should also 
decay more slowly at large r, but this effect is too subtle to be visible with our statistics. 

We next examine the lattice spacing dependence of the correlator at fixed light quark 
mass ratio. Comparing the local correlators C(r) obtained on ensembles at different lattice 
spacing is complicated because sampling is naturally done on a lattice scale. Rather than 
rebinning the data to a common physical scale, we compute the partial integral Xti'f') of 
Eq. (22) and plot it in physical (ro) units in Fig. 5. As r increases from the origin, we see 
a peak at short distance coming from the regulated contact term followed by a decrease 
coming from the negative correlator. The onset and width of the peak is determined by 
the effective radius of the topological charge density operator, which is fixed in lattice units. 
Thus as the lattice spacing decreases, the expected negative singularity in the correlator 
is exposed, and the peak increases in height and decreases in width. 

At large r the data approach the asymptotic value of the full susceptibility. The figure 
shows both the integrated raw data and the integral with the fit values for r > rc replacing 
the raw data. The lower panel enlarges the asymptotic region to show the variance reduction 
achieved by the fit. The result also shows a plausible convergence of the asymptotic value 
in the continuum limit. 

Now we point out a practical issue relevant to future extensions of this work, namely, 
whether the topological susceptibility, defined by integrating the correlator of the regulated 
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spacing 


a (fm) 


ro/a (Jcorr 




coarse 


0.12 


3.82 1.8 X 10" 


-4 


fine 


0.09 


5.40 2.5 X 10" 


-4 


superfine 0.06 


7.73 3.3 X 10" 


-4 


ultrafine 


0.045 


10.39 4.4 X 10" 


-4 



TABLE IV. Error o"corr in Xtifo), the short-distance contribution to the topological susceptibility, 
at sea quark mass rriud = 0.2ms for various lattice spacings. The error is adjusted to the same 
sample size, autocorrelation, and lattice volume. 

topological charge density operator, has a feasibly accessible continuum limit. This will 
be the case if the variance in the integral of the correlator for fixed physical volume and 
statistical sample size does not diverge as the lattice spacing decreases. We examine Xt(^) 
at a fixed physical distance r as the lattice spacing decreases. For r < ro/2 we find that 
the variance actually decreases for a G [0.045,0.12] fm. But for such a small range in r, 
the behavior of the integrated correlator is strongly influenced by the size of the topological 
charge density operator. The larger radius r = ro is safely outside the width of the operator 
and in a region where, for a G [0.045, 0.12] fm, the integrated correlator Xt(^) is well past the 
peak, as we can see from Fig. 5. We show the error in Xt(^o) as a function of lattice spacing 
in Table IV. This statistical error is adjusted for autocorrelations, sample size (factor of 
yiV/iVo), and lattice volume (factor of ^Jv/Vo) for A^o = 500 and I/q = 100 fm^ We see that 
the adjusted error grows approximately as 1/a over this range. This trend suggests that it 
will be increasingly expensive to push to smaller lattice spacing with our scheme. However, 
the continuum limit is nonetheless finite, and our results demonstrate that the method gives 
reasonable errors over the range of lattice spacings considered. 

C. Topological susceptibility 

Our results are summarized in Table VI and Fig. 6. Since chiral perturbation theory 
predicts the behavior as a function of the mass of the taste-singlet pion, we also list estimates 
of that mass. Unlike the Goldstone pion mass, the mass of the taste singlet is not measured 
directly on all of our ensembles. However, to a good approximation, splittings of the squared 
masses of the pion taste multiplet are known to be independent of the light quark mass at 

19 



spacing 


a (fm) 


rgAM2 


coarse 


0.12 


1.136 


fine 


0.09 


0.437 


superfine 


0.06 


0.143 


ultrafine 


0.045 


0.087 



TABLE V. Mass splittings (difference in squared masses) between Goldstone and taste singlet 
pions 

fixed lattice spacing [36]. So if the splitting is measured for one light quark mass for a given 
lattice spacing, the taste-singlet pion mass can be reconstructed from the Goldstone pion 
mass for other light quark masses at the same spacing. Table V lists the splittings for the 
categories of lattice spacings in this study. They were used to obtain the values in Table VI. 
The largest error in the estimated splittings is less than 5%, which bounds the error in the 
abscissa of the plot. We have chosen Tc to lie within the fit range. We have found that 
within this range our results vary by less than one standard deviation. 

D. Continuum extrapolation 

To model a continuum extrapolation, we fit our data to the following form: 

= Co + Ci{a/rof + [c2 + c^ia/r^f + Ci{a/rQf]/{m^jrof . (24) 

This model assumes that lattice artifacts scale as . The fit yields x^/df = 8.8/13. In 
Fig. 6 the resulting fit curves are shown, and three representative points in the continuum 
extrapolation are also plotted. Also plotted is the prediction of Eq. (11) using f.„ro = 
130 X 0.454 MeV-fm with and without our continuum-extrapolated asymptotic quenched 
topological susceptibility Xtf'o ~ 0.0523(29). From the fit itself we obtain f.^ = 132(6) MeV, 
which is better than expected for tree-level chiral perturbation theory. 

IV. CONCLUSIONS 

We have presented an extensive study of the topological susceptibility on 18 (2 + 1)- 
fiavor asqtad lattice ensembles and three quenched lattice ensembles. The susceptibility is 
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10/5' 




range 


(a) 


rja 


'^o'^li iXt<)ro 


iXt>)ro 


Xtro 










coarse 






6.85 


0.05/0.05 


[8.0, 


12] 


10 


4.746 0.0461(14) 


-0.0006(2) 


0.0455(14) 


6.83 


0.04/0.05 


[8.0, 


12] 


10 


3 997 0422(13') 


-0.0008(2) 


0.0414(13) 


6.79 


0.02/0.05 


[8.0, 


12] 


10 


2.580 0.0364(10) 


-0.0009(1) 


0.0355(10) 


6.76 


0.01/0.05 


[8.0, 


12] 


10 


1.872 0.0315(08) 


-0.0015(1) 


0.0300(08) 


6.76 


0.007/0.05 


[8.0, 


12] 


10 


1.665 0.0309(09) 


-0.0015(1) 


0.0294(09) 


6.76 


0.005/0.05 


[8.0, 


12] 


10 


1.517 0.0289(07) 


-0.0021(1) 


0.0267(07) 


8.00 


quenched 


[6.0, 


10] 


10 


- 0.0733(08) 


0.0000(0) 


0.0598(10) 












fine 






7.18 


0.031/0.031 


[10.0, 


18] 


13 


3.626 0.0321(13) 


-0.0018(6) 


0.0303(15) 


7.11 


0.0124/0.031 


[10.0, 


18] 


13 


1 688 0247(09") 


-0.0017(4) 


0.0230(11) 


7.09 


0.0062/0.031 


[10.0, 


18] 


13 


1.074 0.0206(09) 


-0.0031(4) 


0.0174(09) 


7.085 


0.00465/0.031 


[10.0, 


18] 

J 


13 


0.918 0.0188(06) 


-0.0038(2) 


0.0150(06) 


7.08 


0.0031/0.031 


[11.0, 


19] 


13 


0.760 0.0170(06) 


-0.0044(4) 


0.0127(06) 


7.075 


0.00155/0.031 


[12.0, 


18] 


13 


0.601 0.0166(02) 


-0.0047(2) 


0.0118(04) 


8.40 


quenched 


[8.0, 


12] 


10 


- 0.0722(07) 


-0.0000(0) 


0.0593(10) 










superfine 






7.48 


0.0072/0.018 


[12.0, 


25] 


20 


1.177 0.0167(09) 


-0.0023(2) 


0.0144(09) 


7.475 


0.0054/0.018 


[12.5, 


25] 


20 


0.920 0.0148(09) 


-0.0025(2) 


0.0123(09) 


7.47 


0.0036/0.018 


[12.5, 

L ■ 


25] 

J 


20 


0.666 0.0113(09) 


-0.0032(2) 


0.0081(09) 


7.465 


0.0025/0.018 


[13.0, 


25] 


20 


0.510 0.0107(07) 


-0.0037(2) 


0.0070(07) 


7 46 


0018/0 018 


[13.0, 


25] 


20 


0.408 0.0100(05) 


-0.0040(2) 


0.0060(05) 


8.80 


quenched 


[15.0, 


21] 


15 


- 0.0680(06) 


-0.0001(2) 


0.0561(12) 










ultrafine 






7.81 


0.0028/0.014 


[16.0,32] 


27 


0.634 0.0111(10) 


-0.0030(1) 0.0080(10) 



TABLE VI. Fit ranges and cut radius in lattice units and results for the topological susceptibility. 
Also shown are the computed or estimated taste-singlet squared pion masses in tq units and the 
contributions to the total topological susceptibility for distances less than {xt<) and greater {xt>) 
than the cut radius. 
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FIG. 4. Topological charge density correlator vs. r in units of tq for a set of fine lattice ensembles 
(a ~ 0.09 fm) with varying light sea quark masses rriud- Upper: overview. Lower: detail. 
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FIG. 5. Upper panel: integrated topological density correlator Xf(^)'^o ^-^^ '^I'^Q fixed light 
quark mass rriud = 0.2ms for the lattice spacings indicated. Lower panel: detail of the asymptotic 
behavior. The full topological susceptibility is the value at the largest r. The plotted points 
give the result from the raw data without variance reduction. Errors include the adjustment for 
autocorrelations listed in Table III. The solid black curves show the central value of the integrated 
contribution with the fit values replacing the raw data for r > r^- (Values of rc and fit ranges are 
given in Table VI.) The fit curves for a = 0.06 and 0.045 fm are, accidentally, nearly coincident. 
Statistical errors on the solid lines are shown on the right edge of the right panel. They have also 
been corrected for autocorrelations. The fit error for the smallest lattice spacing has the largest 
error bar. The improvement in variance is evident. 



23 
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S-H 

X 

0.02 
0.00 

1 2 3 4 5 

FIG. 6. Topological susceptibility vs. the squared taste-singlet pion mass in units of the Sommer 
parameter tq w 0.454 fm [32]. The brown curve labeled "L.O. 2 + 1 + mo shows the prediction of 
tree-level continuum chiral perturbation theory from Eq. (11) with = 130 MeV, and the dashed 
brown line labeled "L.O. 2 + 1" shows the same prediction without the last term of Eq. (11). The 
remaining curves are fits to the model of Eq. (24). The solid black line is the central value of the 
continuum extrapolation of that model and three representative points on the curve indicate the 
one sigma error. 

defined as the integral of the correlator of the topological charge density. The topologi- 
cal charge density is constructed from a discretized version of FF with smearing to help 
regulate ultraviolet fluctuations. To reduce the variance from large distances, we replace 
the measured values of the correlator at large r by a fit model that builds in the expected 
spectral contribution. 

Our method for determining the topological susceptibility through an integral of the 
topological charge density correlator avoids singularities at zero separation by smearing the 
charge density operator over a fixed local set of lattice sites. A study of the variance in 
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the small- distance contribution suggests that as the lattice spacing is decreased the variance 
grows. At our level of statistics and for the range of lattice spacings we consider in this 
study, this growth is manageable. 

Over the range of lattice spacings and masses in this study, within statistical errors, 
we find good agreement with tree-level staggered chiral perturbation theory and in the 
continuum limit with tree-level continuum chiral perturbation theory, in both cases with 
the expected number of flavors. This agreement supports the assertion that the fourth-root 
procedure for staggered fermions results in the correct number of sea quark species in the 
continuum limit. 
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Appendix A: Ensembles studied 

We use gauge field ensembles generated by the MILC collaboration [5, 36, 37] using 2-1-1 
flavors of improved (asqtad) staggered sea quarks with various light quark masses. Relevant 
parameters of the gauge fleld ensembles in this study are listed in Table VII. They fall 
into four groups according to the approximate lattice spacing, namely coarse (0.12 fm), flne 
(0.09 fm), superflne (0.06 fm), and ultraflne (0.045 fm). The table shows the inverse lattice 
spacing in units of Sommer parameter tq. The pion and ss pseudoscalar meson masses are 
shown in lattice units. 
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10/ff2 


volume 


rriud/ms 




amgs 


ro/a 










coarse 








6.85 


20^ X 64 


0.05/0.05 


0.48454(19) 


0.48454(19) 


3.921 


364 


6.83 


20^ X 64 


0.04/0.05 


0.43488(21) 


0.48647(22) 


3.889 


340 


6.79 


20^ X 64 


0.02/0.05 


0.31134(17) 


0.49012(18) 


3.860 


469 


6.76 


20^ X 64 


0.01/0.05 


0.22439(20) 


0.49427(18) 


3.822 


644 


6.76 


20^ X 64 


0.007/0.05 


0.18903(17) 


0.49324(16) 


3.847 


435 


6.76 


24^ X 64 


0.005/0.05 


0.15970(12) 


0.49261(14) 


3.865 


317 


8.00 


20^ X 64 


quenched 






3.881 


400 








fine 








7.18 


28^ X 96 


0.031/0.031 


0.32003(18) 


0.32003(18) 


5.580 


447 


7.11 


28^ X 96 


0.0124/0.031 


0.20638(18) 


0.32585(17) 


5.420 


509 


7.09 


28^ X 96 


0.0062/0.031 


0.14777(12) 


0.32698(8) 


5.401 


531 


7.085 


32^ X 96 


0.00465/0.031 


0.12851(12) 


0.3269(2) 


5.399 


1000 


7.08 


40^ X 96 


0.0031/0.031 


0.10538(6) 


0.32744(8) 


5.394 


489 


7.075 


64^ X 96 


0.00155/0.031 


0.0750(2) 


0.3275(1) 


5.398 


890 


8.40 


28^ X 96 


quenched 






5.446 


416 






superfine 








7.48 


48^ X 144 0.0072/0.018 


0.13187(8) 


0.20830(12) 7.722 


601 


7.475 


48^ X 144 0.0054/0.018 


0.11420(9) 


0.2075(1) 


7.722 


618 


7.47 


48^ X 144 0.0036/0.018 


0.09353(6) 


0.20731(6) 


7.732 


611 


7.465 


56^ X 144 0.0025/0.018 


0.07843(8) 


0.20764(8) 


7.726 


518 


7.46 


643 X ;l44 0.0018/0.018 


0.06678(3) 


0.20749(4) 


7.710 


799 


8.80 


48^ X 144 quenched 






7.388 


405 






ultrafine 








7.81 


643 X 192 0.0028/0.014 


0.0712(1) 


0.1583(1) 


10.388 


810 



TABLE VII. Simulation parameters for the lattice ensembles used in this study, including measured 
masses of the Goldstone pion and ss meson, inverse lattice spacing in ro units, and number of 
configurations from the ensemble. For taste singlet pions, see Table V. 
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